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Abstract 

Vortex dynamics in a bilayer thin film superconductor are studied through a Josephson-coupled 
double layer XY model. A renormalization group analysis shows that there are three possible states 
associated with the relative phase of the layers: a free vortex phase, a logarithmically confined 
vort ex- ant i vortex pair phase, and a linearly confined phase. The phases may be distinguished by 
measuring the resistance to counterflow current. For a geometry in which current is injected and 
removed from the two layers at the same edge by an ideal (dissipationless) lead, we argue that the 
three phases yield distinct behaviors: metallic conductivity in the free vortex phase, a power law 
I-V in the logarithmically confined phase, and true dissipationless superconductivity in the linearly 
confined phase. Numerical simulations of a resistively shunted Josephson junction model reveal 
size dependences for the resistance of this system that support these expectations. 

PACS numbers: 74.78.-w, 64.60.-i, 75.10.Hk 

Keywords: KT transition, vortices deconfinement 



1 



I. INTRODUCTION 



Topological defects play an important role in many condensed matter systems. A 
paradigm of this are vortices in systems whose energetics may be described by a single angu- 
lar variable 6(r) that is a function of position r, the simplest example being an XY magnet 
(also known as the planar rotor model 0.) The state of the vortices in two-dimensional 
systems determines important physical properties. For example, at high temperatures free 
vortices in superfluid films and thin film superconductors lead to dissipation. At low tem- 
peratures, the defects are bound into vortex-antivortex (VAV) pairs, yielding a state with 
power-law (quasi-critical) behavior in the correlation function (^e l9 ^e~ l9 ^). While neither 
the bound nor the unbound vortex phase supports long-range order Q], there is a well- 
known thermodynamic phase transition - the Kosterlitz-Thouless transition in which 
the vortices unbind at a critical temperature Tkt- 

The situation is dramatically changed when one introduces an explicit symmetry-breaking 
field that aligns the angles (such as a magnetic field in the XY model). One of us 0, 0, 
investigated this situation recently and found that the vortices have three possible phases: 
a free vortex phase, a logarithmically bound VAV pair phase, and a linearly confined phase. 
Can these phases be distinguished in an experimental system described by this model? This 
is the issue we address in this paper. 

The problem of the resistive transition in a single layer superconductor was studied in 
Ref. It was found that for T > Tkt, a current produces a force on a free vortex, whose 
subsequent motion then induces a voltage along the direction of current. This produces a 
linear I — V curve. For T < Tkt, bound VAV pairs have no net force on them due to a 
current and so create no net voltage drop. A voltage is induced by unbound vortices, but 
these are only induced by the current, leading to a non-linear (power law) I — V curve at 
low currents. In this work, we study the analog of this for a system with explicit symmetry- 
breaking: the voltage response of a Josephson-coupled bilayer thin film superconductor with 
respect to an injected counterflow current. As discussed below, the Josephson coupling acts 
as an explicit symmetry-breaking for the relative phase of the order parameters for the two 
layers, #i(r) — #2( r )- Isolated vortices in the relative phase respond to currents in opposite 
directions in the two layers, and their motion induces an interlayer voltage. Thus, to probe 
vortices in the relative phase, current must run in opposite directions in the two layers; 
i.e., we must have a counterflow current. A geometry in which this is induced results when 
current is injected in one layer and removed from the other at the same edge of the sample, 
which to our knowledge was first studied by Fe rre „ and Prange & in the absence of vortex 
excitations. In this situation, the current penetrates the system up to a characteristic length 
scale - the Josephson length - set by the tunneling matrix element and the superfluid stiffness 
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of the layers. Unlike the case for current in an individual layer, the counterflow current can 
decay because an injected charge carrier in the top layer can tunnel to the bottom layer, 
and exit the system on the same edge of the bilayer from which it entered. 

As we show below, when vortices are introduced, their dynamics generates different be- 
havior for the three phases: metallic conductivity in the free vortex phase, a power law 
I — V in the logarithmically confined phase, and true dissipationless superconductivity (zero 
resistance at finite current) in the linearly confined phase. This last phase is qualitatively 
different from anything that occurs for the voltage response of a single superconducting 
layer. As we demonstrate via simulation below, by measuring the voltage difference between 
the two layers at a single edge, we may distinguish the three phases. 

The organization of this paper is as follows. In Section 2, we introduce our model and 
map it to a Josephson coupled bilayer XY model. In Section 3, the phase diagram is studied 
by duality transformation and a renormalization group (RG) analysis. Section 4 discusses 
the I — V curves and finite system effects. The numerical simulations are presented in 
Section 5, and we conclude with a summary. 



II. THEORETICAL MODEL 

We consider a bilayer thin film superconductor, which can be described by a form of the 
Lawrence-Doniach (LD) model [9.J. The LD free energy may be written as 

F = Jd 2 r E„=i, 2 { A l^l 2 + f l^«l 4 + 2 «lW>n| 2 + hcos(9 x - 2 )}, (1) 

where ipi = \ipi\e iei , ip2 = \ip2\e ie ' 2 are the order parameters for the two layers with co- 
ordinates r = (x,y), and the last term is the Josephson coupling between the two layers. 
Neglecting fluctuations of \ip n \, we obtain a free-energy functional with only the phases 
involved, 

F= J d 2 r |2a(Wi) 2 + 2a(W 2 ) 2 + hcos(9 1 - # 2 )} 

= / <Pr\a[V{9 l + 6 2 )f + a[V(0i - 6 2 )] 2 + hcos{6 1 - 2 )}. ^ 

The free-energy JF does not support vortex excitations, because in neglecting variations in 
the amplitudes of j^l we did not allow for the zeros that are necessarily contained in their 
cores. Two reintroduce these, we replace T with a free energy functional F whose degrees 
of freedom are defined on a (square) lattice. By making the replacement 1 — ^9 2 — ► cos(0), 
we allow for configurations in which plaquettes of the square lattice may contain a non- 
vanishing vorticity. The resulting free-energy, which has a form similar to what expects for 
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a bilayer Josephson junction array, is 



F=<* E<rr>> cos[^(r) - 9 1 (r')} + a £ <rr;> cos[£ 2 (r) - 2 (r')] -/iE r cos[^(r) - 9 2 (r)}, 

(3) 

where < r, r' > refers to nearest neighbor sites. Models of this form have been studied 
previously in [lo| ]. 

The analysis of F is complicated by the fact that, in computing the partition function 
Z = J V6e~ F , there are cosines appearing in the exponential. As is well-known P, lllj ]. 
considerable progress can be made if we replace the Josephson coupling form of the par- 
tition function with a Villain model This essentially involves replacing e Jcos9 with 
E^=_oo e " J{e - 2nm)2/2 wherever it appears in the partition function. The replacement works 
because the effective weighting has the same periodicity as the original Josephson coupling 
form, so that we should retain the same possible phases for the system yj. We obtain 
an important simplification by replacing a cosine form with a Gaussian weighting because 
it ultimately allows one to integrate out the angular degrees of freedom 11]. The cost, 
however, is the introduction of effective integer degrees of freedom, which with some work 
can be understood in terms of the vortex excitations of the system 0, [nj]. Following a 
standard trick for rewriting the partition function as a sum over dual degrees of freedom 
P EH we find the partition function for the Villain form of the free energy may be written 
as Z = E s ,t e- Fv ^ T %.5{V ■ Si(r) - T(r))6(V ■ S 2 (r) + T(r)), with 

Fvm[S) = ^ £ r [S?(r) + S|(r)] + ^ £ r T 2 (r), (4) 

where Si j2 are integer fields lying on the bonds of the lattice in each layer, and T is an 
integer field defined on the lattice sites. 



III. VORTEX PHASE DIAGRAM 

A useful representation 0, 0] of the bond degrees of freedom may be written as S± x = 
dy-rn, Siy = —d x m — A, S 2x = d y n, Si y = —d x n+A, where n, m, and A are all integer fields. 
In terms of these the effective free energy Fd for the partition function Z = J2 m n a e ~ F ° is 




As has been discussed elsewhere, the single layer version of this - in which the middle term of 
E]is essentially absent - may be understood as a model of an interface, with domain walls and 
screw dislocations . In that situation the m field allow us to represent configurations with 
closed domain walls, and the A field introduces open domain wall configurations. For Fd, 
we have closed domain wall configurations separately introduced in each layer by the m, n 
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fields, and the A field, rather than completely eliminating sections of closed domain walls 
as in the single layer case, here allows it to shift between the two layers. Because of their 
close analogy, we will call these "kinks" in the domain walls dislocations in the remainder 
of this paper. Physically, the domain walls may be interpreted as worldline trajectories for 
Cooper pairs, and the dislocations represent interlayer tunneling events. 

One advantage of this representation is that the coupling constants are the inverse of 
those in the original functional. This strong-weak duality makes it convenient to study the 
physics with strong inter-layer coupling. We can obtain further insight by looking at the 
problem in the vortex representation. To do this, we apply the Poisson resummation rule 
on the integer fields m, n, and arrive another representation of the partition function 
in terms of integer fields M, N and A. The energy functional associated with these degrees 
of freedom is 



F- V 2tt 2 K\M + N\ 2 v 2t£KW y ( 1 , bin WAfnW 



-^E ? j^(q)#*(q), 



(6) 



where S is the number of lattice sites, M and iV are the vortex numbers for the two layers, 
H = M - N, K = a/2, Q = {Q x ,Q y ), and Q x = 1 - e~ iq * a ,Q y = 1 - e~ iq v a , with a the 
lattice constant. Physically, we understand the M + N integer field as the vortices of the 
symmetric combination of the original layer phases, whereas H represents vortices for the 
antisymmetric combination. Comparing the terms involving H with Eq. 12 in Ref. 0, we 
see that the energetics of anti-symmetric vortices are identical to those of the single layer 
XY model with external magnetic field. It is thus not surprising that the phases of this 
system are analogous to those found in the latter problem. 

It is useful to note that in the representation of Eq. the partition function shows a 
near duality. This can be seen more clearly if one defines J(q) = Q y A(q). In real space J 
is also an integer field, and if one eliminates A in favor of J in Eq. El then a near symmetry 
is apparent under interchange of J and H and K — > l/4nK in Fy. To fully exploit this 
symmetry, we add a term of the form E c \H(q)\ 2 , after which the duality (K — ► 1/AttK, 
E c <-> 1/h) becomes exact 0|. The integer fields H and J may be respectively interpreted as 
the vortex field and a dislocation field, with the added E c term representing a core energy for 
vortices. The duality in this representation shows that if we can determine what happens 
to the dislocations in one part of the phase diagram, we will know what happens to the 
vortices in another. 

To proceed with finding the phases of this system, we need to perform an RG analysis. To 
do this, we wish to work with continuous rather than integer degrees of freedom. Following 
the standard reasoning that the phases depend on the symmetries of the Hamiltonian and 
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not the precise form of the degrees of freedom ^J, we replace the partition function with 
one of the form Z e ff = J T>9 J T>4> J T>ae~ Heff , with 

Heff = J d^Vfc + os|* + ^|V0 2 -ax\* + (7) 

(-y^ cos27T0! + y< t>1 ) + (-y^ cos2tt0 2 + y (j)2 ) + {-y a cos2vra + y a ). 

Here the last three terms are added to favor configurations in which 0i j2 and a are integers. 
The parameters y^ 2 have the physical interpretation of fugacities for vortices in the indi- 
vidual layers, and for small values they favor configurations with vorticity (0, ±1) in each 
plaquette. The parameter y a can be understood as a core energy for a kink of a domain 
wall into one layer from the other, and as has been discussed elsewhere j^| should be taken 
proportional to 1/h to match the energetics of such kinks in the domain wall model Fjj. 

To identify the phases of the vortex system, we wish to find the fixed points to which 
H e ff flows under the RG. We write the effective Hamiltonian as H e ff = H Q + Hi, with the 
unperturbed Hamiltonian 

H = JdV^IWi + «x| 2 + ^|V0 2 - ax| 2 + ^(|g) 2 + \pa> (8) 

with initial value p(£ = 0) = 27n/ 2 (£), and 

H x = {- y4n cos 27T0! + y^) + {- V<h cos2tt0 2 + y^) + YT=2 ^(-l) n (2na) 2n (9) 

with the initial values y 2n (^ = 0) = y a . (Here e e represents the length scale to which the 
degrees of freedom have been integrated out in the RG.) We perform the RG perturbatively 
in Hi, since y a ~ l//t, this amounts to looking at the RG flows in the region of a strong 
coupling fixed point. 

To lowest order in y^ x , y^ 2 , and y a , the RG flow equations can be shown to take the form 



= -(2n - 2} m „ - 27, 2 L(p,(}m,, + 2 <10) 



where (for p « 1) 

L — *S („) 

with A = it /a and £ = ^jKjh. We note that since p is small throughout the flow range, 
the vortex fugacities y^ and yg are strongly irrelevant, indicating that vortices are always 
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bound in both layers for the parameter regime where our calculations apply, large vortex 
core energies and small 1/h (i.e., large symmetry-breaking field.) This may seem at first 
surprising, since we expect that "symmetric" vortices (i.e., the combination M + N in vortex 
free energy, Eq. |HJ) should undergo a Kosterlitz-Thouless transition. This is true, but in 
our formulation it can be found only at order y fay fa because this involves creating vortices 
in both layers at a given position. In the domain wall representation, this corresponds to 
an operator of the form yfa+fa J d 2 r cos [0i(r) + 02 ( r )] becoming relevant, so that closed 
"double" domain walls - closed domain walls lying in both layers at the same time - cannot 
proliferate for the parameters corresponding to the unbound symmetric vortex state. At our 
order in perturbation theory this has no effect, nor do we expect it to qualitatively: since 
"single" domain walls (i.e., closed domain walls in a single layer, generated by either <\>\ or 
(f>2) are proliferated, the presence or absence of these higher order domain walls should not 
affect the qualitative physics of single vortex unbinding. 

The "mass" parameter p plays a crucial role in determining the meaning of the fixed point 
Hamiltonian. If p — > then the field a may effectively remove arbitrarily large sections of 
the closed domain wall sections generated by fluctuations in the fields <\>\ and 02- This 
corresponds to an unbound phase of dislocations. When p > 0, it becomes prohibitively 
expensive energetically to remove large domain wall segments, and the dislocations remain 
bound. Thus, the initial parameters that divide trajectories for which p(£ — > oo) > from 
those in which p(i — > oo) = represent a boundary separating states in which dislocations 
are paired or deconfined 0,0]. The flow equations [TU] have been studied in detail in Refs. Q 
and U and they demonstrate that a transition from p = to p > at the fixed point indeed 
occurs for a critical value of y a ; i.e., a critical value of h. Thus, for arbitrarily small vortex 
fugacity, there is a dislocation unbinding transition at a critical value of the symmetry- 
breaking field, h c , with the unbound state on the small h side of this, and the bound state 
on the large h side. 

Because of the duality between vortices and dislocations, this implies there must also be 
a vortex unbinding transition. The duality suggests that the unbound vortex phase should 
occur for small h and small E c , as one might intuitively expect. There are then three 
possible phases for the system, an unbound vortex phase and two bound phases, the latter 
being distinguished by whether dislocations are in a bound or unbound state. As has been 
discussed in Ref. these two phases may be described in terms of the vortex degrees of 
freedom by whether the VAV pairs are linearly confined, or logarithmically bound. The 
linearly confined phase may be understood as the natural situation for very dilute vortices. 
Returning to the original XY model with a symmetry-breaking field, a state with a single 
high separated VAV pair has its energy minimized by creating a string of overturned angle 
with finite width of order £ connecting the two topological defects jl4j. (For separations 
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smaller than this, the interaction energy is approximately logarithmic, as in the absence 
of the symmetry-breaking field.) This string arises because the symmetry-breaking field 
introduces a large energy cost for states in which the phase deviates from 8 = over a large 
area. Since a vortex requires a 2tt angular rotation for any path surrounding its core, one 
minimizes the loss in interaction energy with the field by confining the rotation to a narrow 
region. These strings may be understood as a degree of freedom dual to the domain walls 
that arose in the dislocation representation. 

The interpretation of the p > phase comes in part from noticing that is irrelevant, so 
that domain walls are proliferated by thermal fluctuations. Because of the duality, we expect 
the same to be true for the strings in the vortex representation of this state. Physically, 
this means the entropy of the strings overwhelms their energy, so that strings of arbitrarily 
large size may be found in a typical configuration. However, the logarithmic attraction 
between VAV pairs is not screened by the proliferating strings, so they remain bound within 
the intermediate phase. Vortex unbinding only occurs at higher temperatures and/or lower 
E c and/or lower h when the entropy overcomes the logarithmic interaction. A detailed 
discussion of how the RG analysis leads to this picture in the single layer case may be found 
in Ref. 0. 

Finally, we note that the existence of all three phases has been demonstrated in the XY 
model with a magnetic field using a simulation method that directly probes the fluctuations 
of the vortices. 



IV. VORTEX EFFECTS IN COUNTERFLOW RESISTANCE 

In this section, we discuss how the different vortex phases could be distinguished in a 
transport experiment. Because the system will contain "double" vortices (pairs of vortices 
stacked across the layers) which behave in a fashion similar to vortices in a single layer 
system, deconfinement of single vortices would be difficult to see in a measurement the 
I — V curve of the bilayer film as a whole. The problem of competing effects due to double 
vortices can be overcome if we consider a counterflow current, which creates forces on a single 
vortex but not on a double vortex. We thus consider a geometry in which the current J is 
injected and removed in the ^-direction from the two layers at the same edge (see Fig. 1). 
The current distribution in such an experiment has been analyzed in the absence of vortices 
and other fluctuations 8], where it was found that the current within each layer j x and the 
interlayer (tunneling) current j z behave as j x ~ sech(x / Xj) and j z ~ sech(x / Xj)tanh(x / Xj) 
with Xj = V ' K/h the Josephson length. These currents are non-vanishing within a strip of 
width ~ A j, so that vortices entering this region will be subject to forces due to the current. 
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To simplify the analysis below, we will model this by a uniform tunneling current within Aj 
of the edge, and take the currents to be zero deeper inside the sample. 

To understand qualitatively the impact of the vortices, consider one vortex in the relative 
phase field $ = Q\ — d%. This configuration can be realized as (M,N) = (1,0) or (0,-1), 
where (M,N) is the vortex number in each layer (see Eq. |UJ). For (M,N) = (1,0), there is 
one vortex in the upper layer, and none in the lower layer. Suppose the current is directed 
in the lower layer so that the resulting motion of the vortex generates a voltage +V relative 
to that in the middle of the system, where we choose the electric potential to be zero. For 
a single vortex of the form (M,N) = (0, —1), the vortex will move in the same direction 
as in the previous case because both the direction of current in the relevant layer and the 
vorticity have changed sign. Since this creates a vortex current of the opposite sign of that 
found in the upper layer, the voltage generated in that layer will have the same magnitude 
on average but opposite sign, —V. Thus, any vortex current induced by the counterflow 
current produces an interlayer voltage at the edge. 

To analyze this in detail, we adopt the approach of Ambegaokar et al. described in Ref. 



15l . which we hereafter refer to as AHNS. AHNS demonstrated that the dynamics of vortices 



in response to a driving current, and the voltage they induce, may be described using a 
Fokker-Planck equation. In this approach, the separation r of a VAV pair is described by 
the stochastic differential equation 

% =-M +»w. < i2 > 

where U is the effective potential for a VAV pair, D is the diffusion constant, and r\ (t) 
represents the thermal noise, with the correlation function < i](t)i](t') >= AD5(t — t'). Since 
the vortices of interest are located in a strip of width Aj near the system edge, we will in 
fact solve a one dimensional problem. From the above stochastic differential equation, one 
obtains a Fokker-Plank equation, describing the relation between T(y), the number density 
of VAV pairs near the edge with vertical separation y, and the vortex current density /. One 
may start from 

F(y, t + At) = j dy'P(y, t + At\y', t)T(y', t), (13) 

where P(y,t + At\y',t) =< 5(y — y(t + At) > y > t is the probability for a VAV pair having 
separation y at time t + At, given that it had separation y' at time t. From Eq. (12) one 
has 

y(t + At) =y>- ^f^At + J t t+M rj^dh (14) 
Expanding P(y, t + At\y', t) to first order in At, we obtain 



P(y,t + At\y',t) 



S(y - y') 

(15) 
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Making use of the fact that f* + t dti f* +At dt 2 < ^(^1)^7(^2) >= 4DAt and integrating by 
parts in Eq. El one can obtain 



which may be rewritten as 

^ = -<§'■ <") 

with 

/ = -2L>exp(^)^[r(y)exp(^)]. (18) 

For a steady-state solution, J is a constant and has the interpretation of the number of 
separating pairs per unit time; i.e. the generation rate for free vortices. The concentration 
of unbound vortices n/ is determined by the balance equation 

In this equation, the first term represents a generation rate for unbound VAV pairs, while 
the second term is due to vortex recombination. In steady state (-^J- = 0), n f ~ 
Finally, we note that the interlayer voltage is generated by a net vortex current, and so is 
proportional to the density of free vortices. Thus we have p ~ rif ~ \fl. Our task is to 
compute /. From Eq. El we have 



r(y )e^) - T(L)e u ^ 
pX ■ 

' dye u ' K - T 



(20) 



' yo 

where L is the system size and yo is a short distance cutoff. At low temperature, to form 
a VAV with large separation, there is a large energy barrier to overcome. Thus to leading 
order in 1/L, we expect T(L — > 00) = 0. So we have 



yo 



dye u ^ 



(21) 



Because there are three phases in this system, there are three effective potentials and result- 
ing resistances to consider. 



A. Linearly Confined Phase 

For the weakest fluctuations (lowest temperatures), we may ignore any screening effects 
of the vortices on their mutual interactions. As discussed above, this leads to a phase in 
which a VAV pair is effectively connected by a string of overturned phase, generating an 
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interaction for highly separated VAV pair that grows linearly with separation. An effective 
potential that describes this situation in the presence of a current j is 



U = e ln(rf a) + e x r j a - jjr 

(22) 

= Eoln(r/a) — Ar, 

where A — jj— £i/a, S\ is the string free energy per unit length, and j is the current. The 
7jr term is introduced to account for the force on a vortex due to the current. Because of 
this, there is a critical current j c = ^y. There are two cases to consider: 

(i) j > jc, A > 0. In this situation there is a critical VAV separation r c = e a/A, for 
which the potential U (r) increases with increasing r for r < r c , and decreases with increasing 
r for r > r c . VAV pairs with separation shorter than r c tend to shrink, while VAV pairs with 
separation larger than r c tend to increase the separation, eventually becoming free vortices. 
For system size L » r c , a saddle point approximation leads to 

T rsj p -U(r c )/k B T 

(23) 

From Ohm's law, the interlayer voltage generated then takes the form 

V = pj~j(j-j c y°/ 2kBT . (24) 

(ii) j < jc, A < 0. In this regime, the potential U(r) increases monotonically to infinity. 
Thus for infinite system, the VAV pairs have to overcome an infinite high barrier to become 
free vortices. We expect no free vortices for an infinite system. In finite size systems, we 
expect corrections to this, which may be estimated as follows. 

I —L 

r /kgT A/kgT l [ rA /k B T ( ^(s /k B T-l), 

-A [ a } 1x0 L V 

- T^-l 1 (25) 

K BJ- ^y Q /k B T c -AL/k B T 

—A a 

^ -A / a\ £n /fc B T c AL/fc B T 

k B T y L> 

and 



V = jp~ J^ 3c k (j-yo/ 2k BT e -y(jc-j)L/2k B T^ (26) 

When system size goes to infinity, the interlayer voltage disappears throughout the region 
j < j c , and we have true dissipationless superconductivity. This is in contrast to single 
layer superconductors, which generically have a power law I — V . For current very small 
~ 0, the resistivity decreases with increasing of system according to exponential law, i.e. 
p ~ e ~eiL/2ak B T _ ^y e w jjj demonstrate this behavior in our simulations below. 
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B. Logarithmically confined phase 



In this phase, the potential for a VAV pair has the form 

U(r) = £oln(r/a) — jjr. (27) 

There is again a critical VAV separation r c = ^ above which they are free. For very large 
system size L » r c , we can apply a saddle point approximation to the integral in Eq. |2U 
and we find 

V ~ pj ~ ji-W2*flT (28) 

This is exactly the result of AHNS. For fixed finite system size L, we consider the situation 
of very small current j ; — >• to obtain an Ohmic dissipation, 

_i 

r *Y°' kBT (29) 



yo 

a\-(i+e Q /k B T) 



(s /k B T + i)(f-y 



Thus the system size dependence of the resistivity follows a power law with system size, 

p rsj £-(l/2+£o/2fc S T)_ 



C. Free vortex phase 

For the free vortex phase, in the limit of infinite system size, the density of free vortices 
is a constant, and we get a linear I — V; i.e., Ohm's law. For finite system size, some bound 
VAV pairs with separation of order L also contribute to the voltage. Thus with increasing 
L, the resistivity decreases (due to the decreasing number of effectively unbound VAV pairs) 
and saturates to a non-zero value as L — > oo. 



V. NUMERICAL SIMULATION 



To test the results discussed above, we studied the dissipation due to vortices in Langevin 
dynamics simulations. The simulation models a Josephson coupled bilayer superconductor, 
in which we directly simulate only the relative phase degree of freedom. This is an important 
simplification because it cuts in half the computer time needed to collect data for a given 
set of parameters, which is often considerable. The dynamics of the angles in this system 
are simulated by integrating the classical equations of motion, so that our system may 
be understood as a variation of the resistively shunted Josephson junction model [l^ j. As 
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described below, we introduce explicit, ideal (dissipationless) leads, which are weakly coupled 
to each of the two edges of the system at x — ±L x /2. In the y-direction we adopt periodic 
boundary conditions. In our Langevin dynamics simulation, the equations of motion are 
taken to be 

Here <3>(r) include the phase difference for the two layers and the two leads. Following Ref. 



19 . we apply a "busbar geometry" in which there is a single phase variable in each of the 



leads. This eliminates the possibility of vortices in the leads themselves, and allows us to 
focus on dissipation due only to the system. In practice, this might be accomplished by 
using superconducting films for the leads that are much thicker than those of the system, so 
that they will have a much higher superfluid stiffness than inside the system itself. An 
important aspect of the simulation is that the coupling between the system and the leads 
must be weak to avoid disturbing the behavior of vortices that approach the edge. For the 
results reported here, we took Kl = Kr = 0.05 (see below) whereas the coupling within the 
system was taken as K = 1, setting our unit of energy. The effective moment of inertial for 
the spins were also taken to be 1. The viscosity rj was taken to be 0.143 in our simulation. 
£(r) is a random torque satisfying < ((r,t)((r',f) >= 2r]T5 r y5(t — t') with T being the 
temperature of the system (chosen to be 1.2 for the results discussed below.) For a system 
of size N x , N y , our Hamiltonian Hxy takes the form 

H XY = -KJ2 <r , r >>cosMr) - $(r')] - h £ r aw[$(r)] 

~Kl Ef=i cos($ L - $(ax + ajy)) - K R £ji cos($ R - $(N x ax + ajy)) 

with K — 1. A typical run consists of 3.9 x 10 8 time steps. Each time step is 0.08 (in the unit 
of yf / K). We also eliminate the initial 1.5 x 10 7 steps in these runs for equilibration. We 
repeated runs for each set of parameters with several different seeds, allowing us to estimate 
the statistical error. 

In principle, we can find the resistances and I-V curves by injecting counterflow current, 
which in this model would act as a constant force on the lead phase variables $l and 
The induced interlayer voltage is then just given by the Josephson relation, Vl,r oc d^i^/dt. 
However, for low currents such simulations become very challenging because there are long 
time scales involved. To simplify our calculation, instead of calculating the I — V curves 
directly, we calculate the resistance for small current (j — > 0) via the Einstein relation. 
The response of a lead variable $l,r to a DC driving force (current) is proportional to its 
diffusion constant, which can be calculated by the correlation function 

D = limr^ 1 < ($'(* + T) - $'(t)) 2 > t , (32) 

where, for example, $' = $ L — $ , and $ is a s pi n angle at some point deep inside the 
system. Here we choose it to be the middle point of the sample. < ... > t refers to an 
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average over different starting times, which is essentially an average over initial conditions. 
Because in our model current acts as a driving force and voltage is an average velocity of 
the lead phases in response, the diffusion constant for the lead variables are proportional to 
the resistance of the system. 

Our simulation results are summarized in Fig. 2. We show the system size (N y ) de- 
pendence of the diffusion constant (i.e. resistance) for /i=0.0286, 0.143,1.43. Our system 
size is N x = 29 (fixed) and A^=40,50,60,75,85. It is clear that there are three different 
behaviors. The curve with solid dots is for h = 0.0287, which shows a straight line with 
slope very close to -1. Here the dashed line is a reference line with slope exactly -1. The 
resistance R ~ Ny 1 comes purely from the geometry of the sample and implies a constant 
resistivity (i.e., R = pN x /N y ). Thus this result is consistent with the system being in the 
free vortex phase. The curve with solid squares for h = 0.143 is a straight line with slope > 1 
on the log- log plot of Fig. 2, indicating that the system size dependence of the resistivity 
has a power law behavior. This is consistent with our expectations for the logarithmically 
confined vortex phase. Finally, the curve with solid triangle for h = 1.43 shows a marked 
downward curvature on the log- log plot. In the insert, this curve is shown in a log-linear 
plot, demonstrating a nice exponential law for the system size dependence. This is consistent 
with our expectations for the linearly confined phase. Thus our numerical results show that 
for a fixed temperature T > Tkt, with increasing h (interlayer tunneling), the resistance 
has three different possible qualitative behaviors. 



VI. CONCLUSION 

By a combination of analytical analysis and numerical simulation, we have shown that a 
bilayer thin film superconductor supports three vortex phases for the antisymmetric combi- 
nation of the layer phase variables: a free vortex phase, a logarithmically confined phase, and 
a linear confined phase. We argued from a Fokker-Planck analysis that the corresponding 
interlayer I — V curves should show measurably different behaviors: metallic conductivity 
for the unbound phase, a power law I — V in the logarithmically bound phase, and true 
dissipationless superconductivity (for an infinite system) in the linearly confined phase. We 
demonstrated that this behavior is consistent with what is found in numerical simulations. 
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FIG. 1: Bilayer system with current being injected and removed from the two layers at the same 
edge by an ideal lead. 
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FIG. 2: The diffusion constant vs. system size N y in a log-log plot with error bars. N x = 29 is fixed. 
The curve with solid dots is for h = 0.0286, corresponding to the free vortex phase, showing Z)(~ 
R) ~ Ny 1 . The dashed line is the reference line with slope -1. The curve with solid squares is for 
h = 0.143, corresponding to logarithmically confined phase, showing D(~ R) ~ Ny~ a , a = 1.74 > 1, 
a power law N y dependence. The curve with solid triangles is for h = 1.43, corresponding to linear 
confined phase. The insert shows the data for h = 1.43 in a log-linear plot, clearly demonstrating 
the exponential system size dependence. 
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